function fig1(filenumber)
count=filenumber;

    N=readdata(count,1);
    S=readdata(count,2);
    T=readdata(count,3);
%     multigridlevels=readdata(count,4);
%     bound=readdata(count,5);
    time=readdata(count,6);
%     dtnom=readdata(count,7);
    dx=readdata(count,8);
    dz=readdata(count,9);
    mx=readdata(count,10);
    mz=readdata(count,11);
    
    
x=zeros(1,S); z=zeros(1,T);
for s=1:S-1 x(s+1)=x(s)+dx(s); end
for t=1:T-1 z(t+1)=z(t)+dz(t); end
    

    
% mTemp=readdata(count,15);
% [grid,xgrid,zgrid]=bilingrid(mTemp,mx,mz,dx,dz,2,2); ind1=find(grid<1);
% avtemp=0*zgrid;
% for t=1:length(avtemp);
%     avtemp(t)=mean(grid(:,t));
% end
% load colormapfig1a;
% figure(1); clf;
% axes('pos',[0.04 0.54 0.42 0.42]);
% [Zgrid,Xgrid] = meshgrid(zgrid,xgrid); grid(ind1)=nan; pcolor(Xgrid/1e3,(Zgrid-40e3)/1e3,grid); shading flat; set(gca,'ydir','rev'); axis equal; colormap(cmapfig1a); colorbar; axis([0 2000 -40 670]); title(['Temperature (^{\circ}C)']); xlabel('Distance (km)'); ylabel('Depth (km)'); clear('mTemp');
% 
% axes('pos',[0.04 0.05 0.42 0.42]);
% plot(avtemp,(zgrid-40e3)/1e3); hold on; set(gca,'ydir','rev');
% plot(adiabat(zgrid,40e3),(zgrid-40e3)/1e3,'r'); axis([0 1900 0 670]); xlabel('Temperature(^{\circ}C)'); ylabel('Depth (km)');
% legend('Horisontally averaged temperature ','Theoretical isentropic temperature at potential \newline temperature T_p=1315 ^{\circ}C','location','southwest');


load surface_time;
S=length(x);
timehistint=linspace(timehist(1),timehist(end),length(timehist)); dt=timehistint(2)-timehistint(1);

lowcut=0;
surfacehisthf=surfacehist;
qhist=2.5*(temphist(:,:,2)-temphist(:,:,1))/(depthrange(2)-depthrange(1)); meanq=mean(qhist(:));
for s=1:S
    tempsurface=interp1(timehist,surfacehist(s,:),timehistint); tempsurface(1)=tempsurface(2);
    temptemp=interp1(timehist,avtemphist(s,:),timehistint);
    tempq=interp1(timehist,qhist(s,:),timehistint);
    avtemphist(s,:)=bandpass(temptemp,lowcut,inf,dt);
    qhist(s,:)=bandpass(tempq,lowcut,inf,dt);
    surfacehist(s,:)=bandpass(tempsurface,lowcut,inf,dt);
    surfacehisthf(s,:)=bandpass(tempsurface,1/100,inf,dt);
    
    for d=1:length(depthrange)
        temptemp=interp1(timehist,temphist(s,:,d),timehistint);
        temptemp=bandpass(temptemp,lowcut,inf,dt);
        temphist(s,:,d)=temptemp;
    end
end

[TIMEHIST X]=meshgrid(timehistint,x);
axes('pos',[0.54 0.05 0.42 0.42]);
pcolor(X/1e3,TIMEHIST,-surfacehist); shading flat; xlabel('Distance (km)'); ylabel('Time/Ma'); colorbar; title('Variations in vertical surface position (m)')

d=9;
axes('pos',[0.54 0.54 0.42 0.42]);
pcolor(X/1e3,TIMEHIST,temphist(:,:,d)); shading flat; xlabel('Distance (km)'); ylabel('Time/Ma'); colorbar; title(['Temperature variations at ' num2str((depthrange(d)-40e3)/1e3) ' km depth (^{\circ}C)']);
